Simulations on Infinite Size Lattices 
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We introduce a Monte Carlo method, as a modification of existing cluster algorithms, which allows 
simulations directly on systems of infinite size, and for quantum models also at f3 = oo. All two- 
point functions can be obtained, including dynamical information. When the number of iterations 
is increased, correlation functions at larger distances become available. Limits q — > and uj —* 
can be approached directly. As examples we calculate spectra for the d=2 Ising model and for 
Heisenberg quantum spin ladders with 2 and 4 legs. 



Standard Monte Carlo simulations are limited to sys- 
tems of finite size. Physical results for infinite systems 
have to be obtained by finite size scaling, assuming that 
one knows the correct scaling laws, and assuming that 
the data are already in a suitable asymptotic regime, ft 
is therefore very desirable to obtain results also directly 
at infinite system size. We will show how to do so, with 
only a small modification of existing cluster algorithms, 
by using the cluster representation of the models to cal- 
culate two-point functions. In the quantum case we can 
then also simulate directly at f3 = oo and calculate cor- 
relation functions and dynamical greens functions. As 
examples we will show calculations for the classical Ising 
model and for quantum Heisenberg ladder systems. 

The Swendsen Wang cluster method |Q for the classi- 
cal Ising model is based on the Edwards-Sokal-Fortuin- 
Kasteleyn |2||| representation 

e pj{ SiSj -i) = £ p5siSj 6 bi . a + (l-p)6 bij , (1) 

b, j =0.1 

=: J2 Wij(*>*i>&«) (2) 

by =0,1 

for the Boltzmann weight of a spin-pair, with p = 
1 — e~ 2 @ J , which enlargens the phase space of spin vari- 
ables Si by additional bond variables bij. The parti- 
tion function Z = J2{ s } S{6 } W(s, b) with total weight 
W(s,b) = Yliij) Wij{si, Sj, bij) is then simulated effi- 
ciently by switching between representations: given a 
spin-configuration s :— {si}, one generates a bond 
configuration b := {bij} with the conditional proba- 
bility p(b\s) ~ W(s,b), thus creating a configuration of 



clusters of sites connected by bonds bi 



1. Given a 



bond configuration, a new spin configuration is generated 
with probability p(s\b) ~ W(s,b). Because of the factor 
5 SiSj 6b ij ,i in Wij, this amounts to setting all spins of each 
cluster randomly to a common new value, independent 
of other clusters. 

Observables O can be computed either in spin- 
representation as 0(s) or in bond representation as 
0(b) = (T, s O(s)W( S ,b))/(Y, s W(s,b)) (so called 
"improved estimators") H. For two-point functions 



0(a) - 
simple: 



SiSj the bond representation is particularly 



0(b) = <5(sites i and j are in the same cluster) . (3) 

Thus two-point functions and derived quantities (includ- 
ing susceptibility, energy, specific heat) can be computed 
from the properties of individual clusters. 

A variant of the Swendsen- Wang method is Wolff's sin- 
gle cluster method pi]. Given a spin configuration, only 
one of the bond-clusters is constructed, namely the clus- 
ter which contains a randomly chosen intial starting site 
xq. All spins in this cluster are then flipped. One advan- 
tage is that the single cluster will on average be larger 
than the average Swendsen- Wang cluster, so that its flip 
results in a bigger move in phase space. Using eq. (||), 
the correlation function can then be measured as 



C(r) = (si s l+r ) = ( 



- E 

% in cluster 



S(i + r in cluster)) 



(4) 



where V c i is the number of sites in the single cluster, and 
the brackets (...) on the right refer to the Monte Carlo 
average. The susceptibility is then \ = P(V c i), the aver- 
age energy can be measured either as E = —J2dC(l) or 
from the average number of bonds in the cluster 0] , and 
the specific heat can be obtained either as a numerical 
derivative or directly from the bond representation H. 

Our new infinite system method now employs the sin- 
gle cluster method, except that it starts each cluster not 
from a randomly chose site, but always from the same 
site xq (e.g. the origin of the coordinate system). Our 
method satisfies detailed balance and ergodicity (for any 
finite region including the origin) similarly to Wolff's 
method. For the ferromagnetic Ising model, we begin 
with an initial staggered spin configuration of unlimited 
size. (Only a finite part if this configuration will have to 
be stored). We iterate the following two steps: (1) For 
the current spin configuration, a cluster containing site Xq 
is constructed with Swendsen- Wang bond-probabilities, 
and (2) all spins in this cluster are flipped, resulting in a 
new spin configuration. The current cluster is discarded 
after each iteration. 
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After a sufficient total number N(r) of iterations, this 
process will "equilibrate" all spins within a radius r 
around xo , and two-point functions can then be measured 
within this area. With increasing number of iterations, 
the cluster will occasionally reach larger and larger dis- 
tances. Because of eq. (Q), the probability to do so is 
(roughly) proportional to the correlation function C(r). 
A region of the lattice at distance r from xq will become 
"equilibrated" after the cluster has reached that region a 
sufficient number n eq of times, thus N(r) ~ tjjtj- Since 
the cluster can grow without bounds, we are calculat- 
ing correlation functions for the infinite lattice, while we 
only need to store the spin configuration within the area 
actually reached by a cluster during a finite run. 

Fig. |I| shows as an example the correlation function for 
the two dimensional fsing model at j3 — 0.42 < f3 c = 
0.44068.., calculated according to eq. (0). For compari- 
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FIG. 1. Correlation functions of the d = 2 Ising model at 
P = 0.42. 

son we show results for finite size lattices, with periodic 
boundary conditions. Usually, the infinite lattice result 
can only be obtained as the extrapolation of the finite 
lattice ones. In contrast, with our new method we can 
obtain the infinite latttice result in a single simulation, 
producing the lower straight line in fig. EL We checked 
that the method works just as well in three dimensions. 

The calculation of error bars for C(r) needs special 
care. For each distance |r|, we demand that the cluster 
contains pairs of sites at that distance |r| d_1 x n eq times, 
i.e. it reaches each site at that distance 0(n eq ) (= O(10)) 
times, before we consider data at that distance thermal- 
ized. Only then do we start to accumulate measurements 
for C(r) according to eq. (^) for all sites i, j in the cluster. 
Alternative similar procedures are possible. 

Note that the overall computational time is by far dom- 
inated by measurements, both in usual cluster algorithms 
and in our method: In order to obtain a relative error e 
for C{r) , a large number of about I/e 2 clusters with sites 
at distance r need to occur, each of which is generated 
with a probability of about 1/C(r), in both cases. Thus 
our method is at least as efficient as usual cluster al- 
gorithms, while providing results for the infinite lattice 
directly. 

The average computational time per cluster is propor- 
tional to the average cluster size, which is proportional to 



the susceptibility (see below eq. (|j)). The effective sim- 
ulated system size grows with the number of iterations. 

When will the new method work ? The computational 
effort and thus the susceptibility needs to be finite. In 
addition, there must not be a finite probability for a per- 
colating cluster, which implies (3 < (3 C . The present im- 
plementation of the new method will thus work anywhere 
in the unbroken phase, right up to the critical point. Sys- 
tems with long range order may become accessible with 
a suitable modification of the method. We note that in- 
deed points away from the transition are often the region 
of interest, when comparing results of simulations to ex- 
perimental measurements. 

For nonrelativistic quantum systems, the loop algo- 
rithm H provides a cluster method. It is based on an 
enlargened representation in terms of the original spin 
operators and additional loop operators, similar in spirit 
to the Edwards-Sokal-Fortuin-Kasteleyn representation, 
eq. ([j]). We will in the following specialize to the spin i 
quantum Heisenberg model 



H = J 



E 



lists 



+ <?7 



■ srs+) + \s*s; 



(•5) 



The two-point function for the single cluster version of 
the loop algorithm then reads 

(i (S+S- + SrS+)) = (5(i and j on the loop)} 
(SfSj) = (Sf S 7 - S(i and j on the loop)) 

(6) 

Therefore our approach to simulate infinite size systems 
can be used in the same way as for the Ising model. "In- 
finite size" here can be applied to the spatial directions, 
as well as, independently, to the direction of imaginary 
time, yielding simulations directly at = oo, while re- 
taining all dynamical information. Note that our method 
is not a projector method. We can simulate directly at 
L = oo and — oo. This also enables us to obtain 
the limits q — > and u> — > directly from the simula- 
tions. Contributions from q = or lu = 0, which are 
in general different from those limits, and are present in 
normal simulations with periodic boundary conditions, 
are completely avoided in our approach, as well as the 
double limit L — > oo and q — > that is usually required. 
When the number of iterations increases, clusters will 
reach larger distances in spatial and/or imaginary time 
direction, thereby improving the resolution in q and u>. 

Again, our present method will be applicable provided 
that there is no percolation and that the unsubtracted 
correlation function drops off sufficiently quickly (i.e. 
that the corresponding susceptibility is finite). 

As an example we studied spin ladder systems fij] 
with N = 2 and 4 legs for the isotropic antiferromagnet 
(A = 1). We used the discrete time version of the loop 
algorithm (At = 1 /16) . The method can be applied just 
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2. Equal time staggered spatial correlation function 
= n and (3 = oo, for N=2 (top) and N=4 (bottom) 



rA of our data directly provides estimates 
0.5059(4) at N = 2 and A = 0.19(1) at 



as well in continuous time |9| . Fig. |2j shows results for the 
equal time staggered spatial correlation functions along 
the chains, which behave similarly to the Ising case. A fit 
to the infinite lattice result gives £ = 2.93(2) for N = 2 
and £ = 8.2(1) for N = 4. Fig. | shows greens functions 
for L = oo, the infinite size system. Whereas finite tem- 
perature calculations give results periodic in imaginary 
time, which have to be extrapolated, the new approach 
provides the /? = oo (T = 0) result directly. With in- 
creased number of iterations, the infinite system results 
become available at larger distances both in spatial and 
in imaginary time direction. A fit to the exponential de- 
cay G(t) ~ e 
for the gaps A 
N = 4, consistent with previous results ||. We also show 
results for L = 20 and j3 = oo to exemplify the effect of 
finite size systems. 

Continuing the imaginary time greens function to real 
frequencies by the maximum entropy method provides 
the spectra of fig. ^, in which the gaps, the single magnon 
peaks, and higher excitations for N — 4 are clearly visi- 
ble. This appears to be the first time that the spectrum 
for N — 4 has been calculated. 

Let us compare our method to other approaches. Ex- 
act diagonalization provides ground state results but is 
limited to small systems. Projector methods are also 
limited in system size and rely critically on proper con- 
vergence to the ground state. DMRG can achieve = oo 
(mostly without dynamical information) for fairly small 
systems, or infinite size for large temperatures and finite 
At. The most powerful method to extrapolate to infinite 
system size is the Finite Size Scaling method of Kim fic| l 
and Caracciolo et al. JllJ , which allows extrapolation at 
correlation lengths far larger than the system size, but re- 
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FIG. 3. Greens functions {S(q, 0) S(q, t)) at q — (tt,tv) for 
L = oo, the infinite size system, with N = 2 (top) and N = 4 
(bottom). Results at L = 20, [3 = oo have been added to 
exemplify finite size effects. 
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quires knowledge about scaling and corrections to scaling 
of the model. 

In summary, we have introduced a new method, as a 
small modification of existing cluster methods, to simu- 
late both classical and quantum systems at infinite sys- 
tem size and/or at zero temperature, while obtaining all 
two-point functions and derived quantities. Larger dis- 
tances (in space and/or imaginary time) are accessed for 
longer simulations. The method in its present form is 
applicable to all cases where existing single-cluster algo- 
rithms can be used, as long as the unsubtracted two-point 
correlation function decays sufficiently quickly, including 
any (3 < (3 C in the unbroken phase. Such points away 
from a transition are often the region of interest when 
comparing simulation results to experimental measure- 
ments. 
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